# Script name  : Relationship_infer.R
# Author       : Xinwu
# Email        : lexinwu@genenergy.cn
# Describe     : Calculated Phi value and inferring pairwise relationships from plink identity-by-descent report.
# Created  Time: 2024-01-24 09:48:01
# Finished Time: 2024-01-24 16:58:06
# License	     : (C)Copyright 2022-2024, Genergy Biotechnology (Shanghai) Co.,Ltd.

if (!require("argparse", quietly = TRUE))
	install.packages("argparse", repo = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
usage  <- "Rscript Relationship_infer.R --IBD_result relationship_infer_result.genome --out_dir ./"
parser <- ArgumentParser(usage = usage)
parser$add_argument("--IBD_result", action = "store", type = "character", required = TRUE, help = "plink identity-by-descent(IBD) report, file end by .genome")
parser$add_argument("--out_dir",    action = "store", type = "character", required = TRUE, help = "Where the results save.")
args <- parser$parse_args()
options(scipen = 10, digits = 5)
suppressPackageStartupMessages(library(magrittr))

# remove "/" if out directory end of "/"
if(substr(args$out_dir, nchar(args$out_dir), nchar(args$out_dir)) == "/"){
  args$out_dir <- substr(args$out_dir, 1, nchar(args$out_dir) - 1)
}
if(!dir.exists(args$out_dir)){
  dir.create(path = args$out_dir, showWarnings = FALSE, recursive = TRUE)
}
if(!file.exists(args$IBD_result)){
  cat("[ERROR: ] Please check IBD_result is sure exists!\n")
  q()
}

data <- data.table::fread(args$IBD_result, header = TRUE)
data$Phi <- (data$Z1 / 4) + (data$Z2 / 2)
data$Relationship <- "UN"
data$Relationship[data$Phi > 0.354] <- "MZ"
data$Relationship[(data$Phi > 0.177) & (data$Phi < 0.354) & (data$Z0 == 0)] <- "PO"
data$Relationship[(data$Phi > 0.177) & (data$Phi < 0.354) & (data$Z0 != 0)] <- "FS"
data$Relationship[(data$Phi > 0.0884) & (data$Phi < 0.177)] <- "2nd"
data$Relationship[(data$Phi > 0.0442) & (data$Phi < 0.0884)] <- "3rd"

# 在样本名前加 ' 因为样本名是纯数字且比较大的话, excel 打开会变成科学计数法.
data$IID1 <- gsub(pattern = "^", replacement = "'", x = data$IID1)
data$IID2 <- gsub(pattern = "^", replacement = "'", x = data$IID2)

data <- subset(data, data$Relationship %in% c("MZ", "FS", "PO", "2nd", "3rd")) %>% 
	    dplyr::select(IID1,IID2,Z0,Z1,Z2,Phi,PI_HAT,Relationship)
write.table(x = data, file = paste0(args$out_dir, "/Relationship_inference_result.xls"),col.names=T,row.names=F,quote=F,sep = "\t")


# 依据
# an estimated kinship coefficient(Phi) range > 0.354  MZ
# an estimated kinship coefficient(Phi) range [0.177, 0.354]  1st-degree(FS PO)
# an estimated kinship coefficient(Phi) range [0.0884, 0.177] 2nd-degree
# an estimated kinship coefficient(Phi) range [0.0442, 0.0884] 3rd-degree

# 如何区分 FS/PO ? 再引入 Z0
# an estimated kinship coefficient(Phi) range [0.177, 0.354]  and Z0 = 0  PO
# an estimated kinship coefficient(Phi) range [0.177, 0.354]  and Z0 > 0  FS